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ABSTRACT 

We consider the stability of warping modes in Keplerian discs. We find them to be 
parametrically unstable using two lines of attack, one based on three-mode couplings 
and the other on Floquet theory. We confirm the existence of the instability, and in- 
vestigate its non-linear development in three dimensions, via numerical experiment. 
The most rapidly growing non-axisymmetric disturbances are the most nearly axisym- 
metric (low m) ones. Finally, we offer a simple, somewhat speculative model for the 
interaction of the parametric instability with the warp. We apply this model to the 
masing disc in NGC 4258 and show that, provided the warp is not forced too strongly, 
parametric instability can fix the amplitude of the warp. 

Key words: accretion, accretion discs - hydrodynamics - instabilities - turbulence 
- waves. 



1 INTRODUCTION 



Thin discs are an important model for accretion flows in a 
number of difi'erent astrophysical objects, including active 
galactic nuclei, interacting binary systems, and young stars. 
Most disc models are flat, i.e. the disc is planar, for sim- 
plicity and because there has been little theoretical or ob- 
servational motivation to consider non-coplanar, or warped 
discs - although observations of the accreting neutron star 
system Her X-1 have long been interpreted as requiring a 
precessing disc that lies out of the plane of the binary sys- 
tem (Tananbaum et al. 1972; Katz 1973; Roberts 1974). 

A pair of developments in the last few years now 
strongly motivates theoretical interest in warped discs. First, 
Pringle (1996) found an instability of a thin disc model that 
tends to produce large scale warps. The instability is driven 
by radiation torques exerted on the disc as it is differentially 
illuminated by the central source and then re-radiates. Sec- 
ond, the megamaser in the nucleus of NGC 4258 has been 
resolved by the VLBA and can be understood as an almost 
precisely Keplerian, but warped, thin disc. The maser spots 
trace out the position of the disc and allow accurate mea- 
surements of its distance, the mass of the central object, 
and the parameters of the warp (Miyoshi et al. 1995). Of 
course, warps may also be excited by direct forcing, e.g. by 
tidal interaction of a T Tauri disc with a binary companion 
(Papaloizou & Terquem 1995). 



Warps in unmagnetized, non-self-gravitating, inviscid 
Keplerian discs are rather special, as we shall review in §3 
below. This is because a warp (with azimuthal wavenum- 
ber m = 1) that is quasi-static in an inertial frame sets up 
vertically-varying horizontal pressure gradients in the disc 
that oscillate at the rotation frequency f2 as viewed in a lo- 
cally corotating frame. In Keplerian discs this is resonant 
with the epicyclic frequency, and so modest warps can ex- 
cite large amplitude motions - we shall call them in-plane 
oscillations - in the disc. As a result, warps and bending 
waves propagate approximately non-dispersively, at nearly 
the sound speed, in a Keplerian disc. In a non-Keplerian 
disc, the in-plane oscillations are much weaker and bending 
waves are dispersive (Papaloizou & Lin 1995; Ogilvie 1999). 

The special nature of warps in Keplerian discs raises a 
number of interesting astrophysical issues. Are the in-plane 
oscillations stable? We shall show via two different lines of 
analysis that they are not in §3 and §4. For small warp am- 
plitudes, the warp is unstable to a parametric instability, a 
possibility presciently hinted at by Papaloizou & Terquem 
(1995). Next, what is the non- linear outcome of the insta- 
bility? This we explore numerically in §5, which examines 
the immediate decay of an initially unstable warp. Finally 
we discuss the implications of this for the disc of NGC 4258 
in §6, using a simple model for the interaction of the warp 
with the parametric instability. 
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2 BASIC MODEL 

To analyse this putative instability of warps in discs, we need 
to formulate an equilibrium to perturb about. In doing so we 
shall make two seemingly severe approximations which we 
believe do not affect the existence or nature of the instability. 

First, we shall consider a radially local model of the 
disc, the well known shearing-sheet model (Goldreich & 
Lynden-Bell 1965). In this approach the equations of motion 
are expanded to lowest order in H{r)/r {H being the disc 
scale height) about some point corotating with the disc at 
(r, 0, z) — {ro, (jjo + ^t, 0). One then erects a local Cartesian 
coordinate system [x, y, z) = {r — ro,ro{4> — 4>o — ^t), z) ro- 
tating with uniform angular velocity = Q,z. In this limit, 
the momentum equation is 

^ = --VP-2nxv + 2qn'^xx-Q?zz. (1) 

The equilibrium velocity field is then Vy^ = —qQ,x, where 
q = — dln57/dlnr = 3/2 for a Keplerian disc. The continuity 
equation takes its usual form, and we shall, unless otherwise 
stated, assume an isothermal equation of state, P = c^p. 
The equilibrium density state is then p = po exp(— 2^/2//^), 
with H = Cs/n. 

Notice that we have neglected magnetic fields com- 
pletely. In a magnetized disc turbulence will develop due 
to the Balbus-Hawley instability (Balbus & Hawley 1998), 
and the analysis presented here will be invalid. In that case 
one must study the interaction of the resulting MHD turbu- 
lence with the warp, a project which, at present, can only 
be undertaken numerically (Torkelsson et al. 2000) . In discs 
that are weakly magnetized the linear hydrodynamic anal- 
ysis that follows will still be applicable over a limited range 
of wavevectors, although it seems likely that Balbus-Hawley 
initiated MHD turbulence will ultimately dominate. There 
may also be conditions under which astrophysical discs are 
so poorly ionized that they do not couple to the magnetic 
field (cf. Gammie 1996; Gammie & Menou 1998). 

Our second approximation, which applies only to the 
analytical approaches used below and not to the numerical 
experiments, is axisymmetry for the bending waves and for 
the parasitic instabilities that feed on them. This may seem 
peculiar as the most astrophysically interesting warps have 
m = 1; such modes have the unique property that they can 
be stationary in an inertial frame, and furthermore are what 
is observed in the unique megamaser disc in NGC 4258. Nev- 
ertheless m — 1 modes are "almost axisymmetric" in that 
{l/r)d^ ~ 1/r, so it is consistent with the local approxima- 
tion to neglect terms of this order and consider only axisym- 
metric bending waves for the basic state. To the extent that 
such terms are negligible, an m = 1 mode with frequency 
u) — and an axisymmetric mode with uj = Q are locally 
indistinguishable because both have the same frequency in 
a frame corotating with the fluid. That we consider only ax- 
isymmetric parasitic instabilities on this basic state is per- 
haps a more questionable approximation, but one that we 
relax in the non-linear regime. 



& Pringle 1993; see also Goodman 1993; Ogilvie & Lubow 
1999). We assume from the outset that both the vertical 
structure of the disc and its dynamical response are isother- 
mal, so that J — 1 and the Brunt- Vaisala frequency van- 
ishes. 

Consider small, axisymmetric Eulerian perturbations 
{5v,SP) from the basic state, with time and space depen- 
dence exp[i{kx — ijt)] understood. These perturbations obey 



-iM^ 
p 



5P_ 
P 



■ , SP \ ^2 , 
lU) I \ — il Z OVz 

P 



(2) 
(3) 
(4) 



Here k is the epicyclic frequency, given by = (4-2g)r2^ 
One can then obtain a single second-order equation for the 
pressure perturbation. 



di 



t) 



1 - 



p 



(5) 



The eigenfunctions are Hermite polynomial^ (e.g. 
Abramowitz & Stegun 1965), 



5P 



— «He„(C) = (-l)V 



(6) 



P dC" 

where C, = z/H and n is a non-negative integer. We recall 
their differential equation 



He::(C)-CHeUC)+nHe„(C) = 
and orthogonality relation 



(7) 



He^(C)He„(C)dC = n!V2^5„„ 



(8) 



The corresponding eigenfrequency tj„ therefore satisfies 



= n G {0,1,2,3,...}. 



(9) 



For future reference we also note the corresponding La- 
grangian displacement. 



ikH 



AA„He„(C), 



e«,„ = 2kH( ."^^ , )AA„He„(C), 
?3,n = A/'„He^(C) = A/'„nHe„_i(C), 
where the normalization 



, 2 



kHuji 



-1/2 



H 



(10) 



(11) 



3 LINEAR THEORY OF KEPLERIAN WARPS 



is chosen according to 



To set the stage for the analysis that follows, we recount the 
theory of axisymmetric waves in isothermal discs (Lubow 



For 7^1, the solutions can also be expressed in terms of 
Hermite polynomials, but with arguments depending on k and uj. 
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where 

M ■- 



j3 

pa X 



(12) 



(13) 



is the mass of the fluid. The eigenfunction ^„ has n vertical 
nodes in and n — 1 nodes in ^z. The dispersion relation 
(^) has two solutions (both positive) for uJ^ at each n: an 
"acoustic" branch with oj^ > c^k^ + and an "inertial" 
branch with < k^. The bending waves correspond to 
n — 1, and the density wave, for which = and dz^x = 
dz(,y = 0, lies on the acoustic branch at n = 0. 

In the special case of a Keplerian disc, n — Q, the 
bending-wave branches cross over at = with non- 
vanishing group velocities ±Cs/2. This happens because of 
the coincidence of the natural frequencies of horizontal and 
vertical oscillations. A disturbance with long wavelength, 
kH ^ 1, is simultaneously close to the Lindblad and verti- 
cal resonance conditions. The form of the dispersion relation 
for small kH is 



(jj ~ r2 ± -Cafc, 



and the normalized eigenfunction in this limit is 



^ = {±iz,±2z,H). 



(14) 



(15) 



The physical velocity field, with arbitrary normalization, is 



Sv = S (^zcos{0,t), — sin(nf), sm{Qt)J 



where 



dSvx 



dz 



(16) 



(17) 



is a measure of the shear rate in the warp. 

This velocity perturbation consists of two parts. The 
vertical motion is simply a uniform oscillation at the vertical 
frequency (which, in a Keplerian disc, is equal to the orbital 

frequency). In a disc with anm — 1 warpj^ it corresponds to 
the fluid following an inclined circular orbit. This motion can 
be eliminated by redefining the frame of reference to follow 
the inclined orbit, with the result that only the horizontal 
components of (^) remain. The horizontal motion is in fact 
an exact solution to the equations of motion in the shearing 
sheet, so S can take on any value. We shall consider the 
stability of this motion; it is in a sense a local model in 
azimuth and radius for the motion in a large scale, m — 1, 
warp. For the special case of a quasi-static m = 1 warp in 
a Keplerian disc, the disturbance can remain close to the 
Lindblad- vertical resonance over the entire radial extent of 
the disc. 

In the limit that the amplitude of the warp is very large, 
so that 5 3> n, the oscillatory character of the motion is 
not important and the problem reduces to the stability of a 
vertically linear shear flow. This problem was considered by 
Kumar & Coleman (1993), who found it unstable. 



t We recall that the same eigenfunctions apply to non- 
axisymmetric modes of small m (m <C r/H) provided that ui 
is understood as the frequency seen in a locally corotating frame. 



4 STABILITY OF WARPS: FLOQUET 
ANALYSIS 

The stability problem is reduced to its simplest possible 
form if we assume that the fluid is incompressible, rather 
than isothermal, and neglect its stratiflcation. The govern- 
ing equations are the same as before, except that V-u = 
and the vertical gravity term is omitted in the equation of 
motion. 

The basic state is now regarded as inflnite in both radial 
and vertical extent, and has uniform density and pressure. 
We divide the velocity field into three pieces: v — vq + vi + 
Sv. The first piece, vq = — ^^xy, is the differential rotation 
of the disc as it appears in the local model; the second piece 
is the "bending wave" , given by the following exact solution 
to the equations of motion (cf. equation |l6| ): 

Vx,! ~ SzCOs[Qt), 



-Sz sin(S7t), 



(18) 
(19) 



and Vz^i — 0. Since the solution is exact, S can be arbitrarily 
large. 

The third piece of the velocity field, Sv, is the pertur- 
bation on the warp whose stability we wish to consider. In 
general, one might think of linearizing the equation of mo- 
tion and projecting Sv on to a suitable basis of functions. 
Note that one cannot find normal modes because the basic 
flow is itself time-dependent. However, a standard technique 
is available for an incompressible fluid if the strain tensor of 
the basic flow is independent of position (Kelvin 1887; Craik 
& Criminale 1986; Bayly 1986). Here one considers a per- 
turbation in the form of a Fourier mode whose wavevector 
evolves in time according to the strain fleld. Specifically, we 
take all perturbed quantities to be time-dependent ampli- 
tudes multiplied by exp[ife(t) ■ x], where 



— = -[\7{vo + vi) 



k. 



(20) 



constant. 



-Skx cos{Q,t), 



(21) 



For an axisymmetric mode {ky = 0) this gives 

dfcz 
dt 

so that kz = kz,o — {S/Q,)kx sm{flt). In words, then, these 
perturbing waves nod back and forth on top of the oscillat- 
ing in-plane motions set up by the warp. In the absence of 
the warp, they would be inertial oscillations (also called r 
modes) with uj^ — Q^kz/k^. 

The amplitudes of the perturbations obey 



— — 2fl Svy + S cos(^U)Svz = —ikx ( — — 
dt \ p 

^ + ^nsvx - ^Ssm{nt)Svz = 0, 



ifc ( ^ 



ikx Svx + ikz Svz = 0. 



(22) 
(23) 

(24) 
(25) 



Note that this is an exact solution since the non-linear self- 
interaction of the mode is zero. 

Recall that kz is time-dependent. These equations can 
be reduced to a single Floquet equation (a linear equation 
with periodic coefficients) of the form 
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^ {k\:^'''5v,) + fit) (k^kz'/^sv,) = 

where 
fit) 



fc2 



^2 ^2 ^2 



(26) 



(27) 



is a function with period 2n/i}. 

First we assume that the shearing motion is small, in 
that S/fl <^ 1. Expanding / in powers of S/fl, 



2fc^ 



lU^nSsin(f7t) + of^) . (28) 



To this order the equation for the evolution of the pertur- 
bation is the Mathieu equation, which has the normal form 



d^ 

+ ri^(a + 2ecosSli)y = 



(29) 



(see, e.g.. Bender & Orszag 1978). 

The Mathieu equation is the classical model equation 
for parametric instability. If e ^ 1, parametric instability 
sets in for a = r? j'^, n = 0, 1, 2, . . .. Here the Mathieu equa- 
tion describes the evolution of an inertial oscillation with 
initial wavevector fco and frequency < O, so a < 1. Instabil- 
ity is then achieved by tuning fco so that n = 1, or equiv- 
alently, so that the inertial oscillation has frequency ±f2/2. 
The inertial oscillation is destabilized by its interaction with 
the warp, which has frequency f2. Near n = 1 the growth 
rate of the instability is eVL. 

Applying these results to the warp, we find instability 
when kx/kzfi = ±\/3, so the wavevector of the perturbed 
velocities is initially oriented at ±60° to the vertical. Then 



}{t) = 



i±^gsin(«.)l+0^^ 



The growing mode has growth rate 



16 



(30) 



(31) 



i.e. the growth rate is proportional to the shear rate in the 
warp. 

If the shearing motion is not small, the Floquet expo- 
nents and associated growth rate can still be obtained nu- 
merically. Instead of using equation ( p^ ) , which is singular if 
5* is so large that k^ passes through zero, one can integrate 
the coupled equations 



d_ 

At 



ASvz 
iSP 



iSP 
P 



fe2 

Ak^kz f iSP 
P 



A;2 



(32) 



over a single period, starting from two linearly independent 
initial conditions. The result of this calculation is that, as 
5* increases, the growth rate s increases less rapidly than a 
simple linear relation. For S/Q = 0.1, 1 and 10, the growth 
rate (maximized with respect to k^/k^fi) is 99.9%, 90.1% 
and 31.1%, respectively, of the value given in equation (^l|). 

We remark that this calculation could easily be ex- 
tended to non-axisymmetric modes, and viscosity could also 



be allowed for. However, the method is restricted to incom- 
pressible fluids. 

In the case of a non-Keplerian disc (k 7^ n) the anal- 
ysis is almost identical. Provided that S still measures 
the amplitude of the a:-component of the epicyclic motion, 
Vx,i = Szcos{Kt), the Floquet equation and growth rate are 
unchanged except that Q is replaced by k. The main differ- 
ence in the case of a non-Keplerian disc is that the in-plane 
oscillations are more weakly coupled to the warp. 

This simple calculation strongly suggests that warped 
discs are hydrodynamically unstable, although it uses a 
model lacking compressibility and the vertical structure of 
the disc. In the next section we remove these limitations, at 
the price of a lengthier and less transparent analysis. 



5 STABILITY OF WARPS: MODE COUPLING 
ANALYSIS 

Parametric instability can be treated as a special case of 
three-mode coupling in which one of the three modes (mode 
1) has a prescribed motion, and hence can be regarded as 
a time-dependent component of the basic state upon which 
the other two modes (2 and 3) are treated as linear pertur- 
bations.^ Sometimes, as is the case here, the wavelength of 
mode 1 is much longer than that of mode 2 or 3. 

In the absence of dissipation, three-mode couplings are 
most easily calculated from an action principle. 



dt d^xCi$,dti,d4,x,t) ^0. 



(33) 



Here indicates a first-order variation with respect to the 
displacements ^, which are regarded as functions of the spa- 
tial coordinates x and time t. Thus x = x + ^{x,t) is the 
actual position of a fluid element whose equilibrium position 
is X. For clarity, we sometimes indicate the basic state with 
an overbar: e.g., p{x) is the mass density in the basic state. 
In the case of axisymmetric motions of the shearing sheet, 
the Lagrangian density to be used in equation (bd) is 



C 



2^2 



u 



}.(34) 



The frequencies fl and k are constants. For simplicity, we as- 
sume that U, the internal energy per unit mass, has the form 
appropriate to an isothermal gas of uniform sound speed c^: 

2 1 



U 



' In p. 



(35) 



The density p differs from p because of changes in volume: 

= 1 + V-l + {d^i^)(d.£,,) - (94^)(9.Cz).(36) 

Therefore, the internal energy (^) contains terms of all or- 
ders in ^. The quadratic and cubic terms are 

1 



U2 
Us 



j(V-^) + dx^zdz^x - dx^xdz^z 
{V-$){dx^xdz^z-~dx^zdz^x 



(37) 
(38) 



■t The amplitude of mode 1 must still be small, 0(t). Our expan- 
sion (equations [37,38]) assumes that the amplitude of modes 2 
and 3 is 0{5€), with 1 > <5 > e. 
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If the spatial derivatives of ^ are small, higher-order terms 
can be neglected. Then U2 and the other quadratic parts of jC 
determine the linear equations of motion, while Us provides 
the dominant non-linear interactions. 

A general disturbance can be expanded in terms of the 
spatial eigenfunctions of the linear normal modes:H 



(39) 



Absent non- linear interactions, the coefficients {qa{t)} 
evolve independently according to 



(40) 



where cOa is the natural frequency of the a* normal mode. 
Non-linear equations of motion for the qs result from 



dtL{{q^},{qo,}) = 0, 



with 
L = 



E 



2 2 ^ 



d^x p{x) Uz 



qcim^ix) +0{q'). (41) 



The spatial eigenfunctions have the normalization (|l^; 
therefore has units of length and qa is dimensionless. 

Consider a situation in which only three modes 
{gi,q2,?3} are significantly excited. Then the Lagrangian 
reduces to 



1 2 



Tqiq2q'i, 



(42) 



where we have divided out an overall factor of MH , and 



dqidq2dqs 



d:'x p(z) Us 



qc.(t)$.^{x) 



(43) 



The behaviour of a system of thr ee oscillators non- 
line arly coupled as in (|4^) is we ll known ( sagdeev & Galeev 
1969 |; [K|umar fc Goodrnan 199e| ) : 

1. The Hamiltonian corresponding to ( ^2[ ) has a local but not 
a global minimum at the origin, because the highest terms 
are odd in q. When, as here, (k2) results from truncation of 
a physical system with a positive-definite Hamiltonian, then 
quartic and perhaps higher terms must be retained if the 
amplitude or coupling is sufficiently large. 

2. In the opposite limit of small amplitudes or weak coupling, 
where the cubic term is small compared to the quadratic 
ones, non-linearity is important only if approximate reso- 
nance exists among the natural frequencies of the three lin- 



S This assumption would be highly questionable for non- 
axisymmetric perturbations: in that case, the modes are probably 
not complete, or if complete, not discrete. 



earized oscillators. Without loss of generalityp], let lui > 
L02 > W3. Resonance exists if wi ~ (j;2 -ftJa, and the accuracy 
of the resonance is characterized by the smallness of 



-. LUl — LU2 — 1^3- 



(44) 



On the other hand, the strength of the coupling can be char- 
acterized by a dimensionless parameter such as 



(45) 



where E is the total energy of the system (^). 
3. For e <C 1, slow exchange of energy among the three os- 
cillators can occur on time-scales ~ (euj)^^ provided |Acij| < 
I eoji I . There is no general restriction on the amount of energy 
that can be exchanged, except that 

AEs _ AE2 _ AEi 
UJ3 1^2 1^1 
and that the energy in each oscillator, 

Ea 



(46) 



/ -2 2 2 > 

[qa +l^aqa, 



(47) 



must remain positive. The selection rule (^6|) can be ex- 
plained, or at least remembered, by saying that a pair of 
quanta of action I = E/u in oscillators 2 and 3 combine in 
pairs to make a single quantum in oscillator 1. 
4. In the special case that all of the energy is initially in oscil- 
lator 1, the energies of oscillators 2 and 3 grow exponentially 
at growth rate 2s, where s is the growth rate of their ampli- 
tudes. The maximum growth rate is achieved when Alu — 0: 



in 



El 



8LjfL02ijJ3 



E2_ 

2E2 



E-i 
2E3' 



(48) 



Of course growth slows once a significant fraction of the 
total energy resides in oscillators 2 and 3. If resonance is 
not exact, then the growth rate is 

s = ^si,, - (Alj/2)2, (49) 

and there is no growth at all where |Aa;| > 2smax. 
5. Finally, if dissipation is present so that the free decay rate 
of the amplitudes of oscillators 2 and 3 is rj, then formula 
( ^9| ) holds with Smax replaced by Smax — V- 

Finding the parametric growth rate due to a large-scale 
bending mode reduces therefore to the evaluation of the cou- 
pling parameter F given by equation ( ^3| ) in terms of the 
linear eigenfunctions and eigenfrequencies {(C^i'^a)}- ^ 
local approximation, the background density p depends only 
upon z, so the [x, z) dependences are separable. Hereafter 
we omit the overbar from quantities of the basic state except 
where needed for clarity. Thus the horizontal dependence of 
the eigenfunction is exp{ikax) . To obtain a non-zero result 
for the coupling (tea), we must have 



ki + k2 + k3 = 0. 



(50) 



In most cases of interest, mode 1 is a large-scale bending 
mode with fci ^ , and modes 2 and 3 have wavelengths 
< 2nH, so to an adequate approximation, k2 ~ — fcs ^ fci, 
and we shall write k for k2. In the limit fci = 0, assuming 
that K = fl, the bending mode is (cf. §3) 

^ We presume that the linearized oscillators are stable, so that 
every ui'^ > 0, and we define lOc, = | \/ oj"^ | 
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(51) 



where S is the vertical shear ol the mode as defined by equa- 
tion ([|). 

Comparing equations (^l|) and (^), one sees that the 
coupling is dominated by terms involving (9z^a;)i, i.e. by the 
vertical strain due to the bending mode. Retaining only the 
leading order terms, 

r^-^ / d^pWc?((v-02(9.c.)3 + (v.03(a.c;)2)^, 

(52) 

where E = J dz p{z) is the surface density. 

Whereas the eigenfunction ( [s]] ) of the bending wave is 
independent of the vertical structure of the disc in the long- 
wavelength limit, the eigenfunctions of the short- wavelength 
daughter modes 2 and 3 are not. These are described in §0 

We are now in a position to calculate the coupling (B^. 
Note that modes 2 and 3 must have different values of n: 
under the reflection z —z, V-^„ has the parity (— )", 
whereas dx£,z has the parity (— )"~^. Therefore F 7^ only if 
?^2 — n^, is odd. This is probably true for any vertical struc- 
ture symmetric about the mid-plane, if n counts vertical 
nodes. For the isothermal gas, the properties of the Hermite 
polynomials are such that F 7^ requires Ins — 712! = 1. We 
must also remember to change the sign of k between the two 
modes, because = ^2 ~ — fcs. Finally, it turns out that the 
horizontal phase of modes 2& 3 must differ by 90° (a fac- 
tor i in our complex notation) in order to get the maximum 
coupling. Taking «: = f2, we find 



n + 



kHuj 



'1/2 



-I -1/2 



(53) 



Here t<;„ and tJ„_i are implicit functions of k and n through 



equation (|9|) and the resonance condition 
For comparison with results obtained in 



-I- LOn-l ~ n. 

M, we consider 



the limit 




kH 


00, 


n 


00, 







It is easy to see from 
and then we obtain 



Foe = 



3%/3 2 
16 ■ 



that this requires kHj^/^n 1, 



(54) 



Taking into account that the total energy of the har- 
monic oscillation ( ^ ) is twice its mean kinetic energy, we 
have E\ = fi'^lgil^ax/S, where |gi|max = is the dimen- 

sionless semi-amplitude of the bending mode. The growth 
rate (Es) reduces to 



IFSI 



(55) 



In particular, for F — > Fo 
3-\/3<S'/32. This is exactly half the growth rate calculated in 



and L02 = ijJ'i = 5^^, 



The explanation appears to be that each daughter mode 
with n ^ 1 and k ~ \/3n can participate in two parametric 
couplings: one with mode (n — 1, —k) and the other with 
(n + 1,— fc). At sufficiently large n, both of the resonance 
conditions ^n.k + '^n±i,-k ~ ^ can be satisfied with nearly 
equal accuracy. As a result, the growth rate of mode (n, fc) 
is doubled. 

It is, of course, not obvious that the mode-coupling cal- 
culation should give the same answer as the Floquet analysis 
of §4, even after this correction. The gas considered here is 
compressible, whereas in §4 it was assumed incompressible. 
As n ^ 00, however, the modes of interest, which belong 
to the "inertial" branch (§3), are effectively incompressible 
near the mid-plane. At high altitudes (z > //) the inertial 
modes are nearly isobaric, and the local contribution to the 
coupling coefficient (^) varies with height. But when modes 
of many (large) values of n and the same |fc| grow together, 
it can be shown that the relative phasing of the mode ampli- 
tudes demanded by parametric instability results in a total 
wavefunction that is concentrated toward the mid-plane. 



6 NON-LINEAR OUTCOME 

The analyses of §§4 and 5 establish the existence of a lin- 
ear instability of hydrodynamic warps. We would now like 
to understand the non-linear outcome of the instability, and 
its non-axisymmetric development, which was ignored in the 
analytic analysis. To do this, we have integrated the hydro- 
dynamic equations in a series of numerical experiments. 



6.1 Initial conditions, boundary conditions, and 
numerical method 

The numerical model retains the basic assumptions used in 
the linear analysis. The parent warp mode is axisymmetric. 
The equations of motion are those appropriate to the local 
model (see equation 1). The disc is Keplerian, &o q = 3/2. 
The fiuid is assumed strictly isothermal. We do not require 
that the system remain axisymmetric. 
In the initial conditions we set 



v^ = Sz + (5, 

vy = -(3/2)J7a; + 5, 

= + 5, 
p = poexp(-^V(2^/')). 



(56) 
(57) 
(58) 
(59) 



Here \&\ < 10"'^ is a uniformly distributed random vari- 
able with mean zero chosen independently for each variable 
in each zone. These perturbations around the equilibrium 
seed the parametric instability; we would otherwise need to 
wait for the instability to grow from machine roundoff error. 
The vertical shearing motion v^c = Sz is the local manifes- 
tation of a large scale warping motion in a Keplerian disc, 
and we shall loosely refer to it as "the warp" . 

The boundary conditions are as follows. We integrate 
the basic equations (^ in a rectangular box of length in 
the radial direction, Ly in the azimuthal direction, and in 
the vertical direction. The grid is centred on x = y = z = G. 
The radial boundaries (at x = ±Lx/2) are subject to the 
"shearing box" boundary conditions (e.g., Hawley, Gammie, 
& Balbus 1995). These require that 
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fi-L^/2, y, z) = /(L,/2, y - qflL^t, z) 
for / = Wa;, Wz, p, where Svy = 



(60) 



+ qQx. The azimuthal 



boundaries (3/ = ±Ly/2) are periodic. The vertical bound- 
aries are reflecting. 

The numerical model is integrated using a version of the 
ZEUS algorithm (Stone & Norman 1992). ZEUS is an ex- 
plicit, finite-difference, operator-split method. The variables 
lie on a staggered mesh, so that scalar quantities are zone- 
centred, while vector quantities are centred on zone faces. It 
conserves mass and linear momentum to machine roundoff 
error. 

The ZEUS algorithm has been extensively tested (Stone 
& Norman 1992; see, e.g.. Stone & Hawley 1997 for a dis- 
cussion of other applications) . Our implementation has also 
been tested on a number of standard problems. It can, for 
example, reproduce standard linear results such as sound 
wave propagation, and standard non-linear results such as 
the Sod shock tube. 

One relevant test of our shearing box implementation 
is uniform epicyclic motion. This can be initiated by setting 
Vx = const, in the initial conditions. The fluid ought then to 
execute epicyclic oscillations with period 27r/(f2(4 — 2qY^^). 
This test is non-trivial: for a naive implementation of the 
Coriolis and tidal forces the amplitude of the epicycle will 
grow or decay. This is because the time-step determined 
from the Courant condition depends on epicyclic phase. 
Growth or decay of the oscillation can be p revented by 
using "potential velocities" v^^x = x/t'I-M^ and Vp^y = 
•\/vy + v'i 1 4 rather than , Vy in the time-step condition. 

A second, relevant test is the evolution of a linear am- 
plitude bending wave. This test was performed in axisym- 
metry. We set Lx,Lz — (10, 10)-ff, and the numerical reso- 
lution Tlx, riz = 64, 64. We introduced a bending wave in the 
initial conditions with wavelength L^- Linear theory gives 
w = 0. 734051 (notice that this is significantly different from 
n, so pressure gradients play an important role in the mode 
dynamics ; the mode is not a mere sloshing up and down of 
the fiuid in a fixed potential) . The measured mode frequency 
(interval between successive zero crossings of the vertical ve- 
locity at a fixed point in the fluid) is a; = 0.7336f2, which 
differs from the true value by 1 part in 2 x 10^, indicating 
satisfactory performance on the test. 



6.2 Non-linear outcome 

The numerical model has seven important parameters: 
Lx,Ly,Lz,nx,ny,nz, and S. Before studying the effect of 
these parameters we will consider a single, fiducial run in 
detail. 

The flducial run has an initial amplitude S = Q. 
The physical size of the box is Lx,Ly,Lz = {A,1Q,Q)H. 
The reflecting vertical boundaries are therefore 3 scale 
heights away from the mid-plane. The numerical resolu- 
tion is Tlx, fly, = 128,128,192, so the resolution is 32 
zones/scale-height. The run ends aX t — lOOfJ^^. 

As diagnostics, we record the evolution of three integral 
quantities. The "epicyclic energy" is 



(61) 



2 


J 1 1 1 1 1 1 1 1 1 1 1 1 





1 






H/Z 


' \ " .i ' 
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' '.'f id--". 
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1 2 



X/H 



Figure 1. The momentum field p5v on a constant- j/ slice through 
the fiducial numerical model. The velocity due to the parent bend- 
ing wave has been removed. 



mode in the absence of an instability. The vertical kinetic 

energy is 



j3 2 
d X pVz . 



(62) 



Absent initial perturbations, _Ek,z = 0. We expect it to grow 
at twice the parasitic mode amplitude growth rate. As a 
measure of the growth of non-axisymmetric structure, we 
define 



NA, 



where 



.{x) 



Ax (c^(x) -l-s^(x)). 



AyAz cos{2mny/Ly)pvx{x). 



(63) 



(64) 



Eepi would be constant for the infinite wavelength warping 



and likewise for Sm (x) . N A^ will grow only if there is a non- 
axisymmetric counterpart to the parametric instability or a 
tertiary nou-axisyinmetric instability that preys on the par- 
asitic modes in the non-linear regime. There are no known 
local, linear, non-axisymmetric instabilities of a Keplerian 
disc. 

Figure 1 shows the perturbed momentum field pSv at 
t = 20Q^^ on a constant-j/ slice. Here 5v is the velocity 
field with the parent bending wave removed. We expect the 
velocities to be dominated by the resonant inertial modes 
which lie (approximately) at 60° to the vertical, and this is 
what is seen in the figure. 

The evolution of i5opi,-Ek,z, and NAi_2 are shown in 
Figures 2, 3, and 4. From Figure 2 it is evident that the 
instability grows rapidly, removing energy from the warp. 
The instability couples strongly to compressive modes, since 
its characteristic velocity amplitude is SH = Cs{S/Q) = Cs- 
As a result shocks are produced through which part of the 
energy is dissipated. The remainder of the dissipation is via 
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Figure 2. Evolution of tlie epicyclic energy in the fiducial nu- 
merical model. 
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Figure 3. Evolution of the vertical kinetic energy in the fiducial 
numerical model. 

numerical averaging at the grid scale. The epicyclic energy 
is reduced by a factor of 5.5 over the course of the run. 

At early times in Figure 2 a small oscillation in epicyclic 
energy is evident. This is caused by the limited accuracy 
of the time integration of epicyclic motions. As a result 
epicyclic energy varies with epicyclic phase. The oscillation 
is stable, however: the epicyclic energy neither grows nor de- 
cays from cycle to cycle. We have confirmed that the ampli- 
tude of this oscillation decreases as zone size and time-step 
decrease. 

Figure 3 shows the development of Sk,z on a logarith- 
mic scale, along with a line showing the analytic prediction 




tn 

Figure 4. Evolution, in the fiducial model, of NAi_2,3, which 
measure lowest-order non-axisymmetric structure (see text for 
definition). The ordinate units are arbitrary. 



for the growth rate. In the linear regime, the upper envelope 
of i5k,z ~ exp(0.51f2t), while the analytic theory predicts a 
growth rate of SVSQ/S w 0.65f2 for the energy. Several fac- 
tors contribute to lower the growth rate below the analytic 
value. We shall discuss these below. 

Figure 4 shows the development of lowest order (m = 
1, 2, 3) non-ajcisymmetric structure in the numerical model. 
The growth rate of NAi is larger (as 0.86f2) than the grovirth 
rate of the axisymmotric instability. This, together with the 
fact that NAi begins to grow only after the parametric in- 
stability roaches the non-linear regime at t w 45S1^^, sug- 
gests that, in this case, growth of non-axisymmetric struc- 
ture is due to a tertiary non-axisymmetric instability, rather 
than the linear development of a non-axisymmetric counter- 
part of the parametric instability. 

Figure 5 shows a colour coded image of density on slices 
through the model in the x — y plane at 2 = (top), the 
X — z plane at j/ = (middle) and the y — z plane at x = 
(bottom). Black is highest density, blue lowest. The slices 
are taken late in the evolution, at t = 80Q~^ . Density vari- 
ations are visible because compressive waves are strongly 
excited. The extended, sharp features are shocks. The r.m.s. 
variation in surface density is ((5S^)^''^/(E) = 0.04. 



6.3 Parameter study 

We have explored the effect of each of the main model pa- 
rameters on the development of the parametric instability 

and its non-linear outcome. A full list of relevant models, 
with model parameters, is given in Table 1. Run 1 is the 
fiducial model. 
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Table 1. List of all numerical models 



No. 


Lx 


Ly 






Uy 


riz 


S 




1 


4 


16 


6 


128 


128 


192 


1 


fiducial run 


2 


4 


- 


6 


256 


- 


256 


0.5 




3 


4 


- 


8 


340 


- 


340 


0.5 




4 


4 


16 


4 


64 


64 


64 


1 




5 


4 


16 


10 


64 


64 


160 


1 




6 


4 


- 


10 


64 


- 


128 


0.5 


periodic, axisymm. 


7 


4 


- 


10 


128 


- 


256 


0.5 


periodic, axisymm. 


8 


4 


- 


10 


256 


- 


512 


0.5 


periodic, axisymm. 


9 


4 


- 


12 


64 


- 


128 


0.5 


periodic, axisymm. 


10 


4 


- 


12 


128 


- 


256 


0.5 


periodic, axisymm. 


11 


4 


- 


12 


256 


- 


512 


0.5 


periodic, axisymm. 


12 


4 


8 


6 


64 


32 


96 


1 




13 


4 


16 


6 


64 


64 


96 


1 




14 


4 


32 


6 


64 


128 


96 


1 




15 


4 


64 


6 


64 


128 


96 


1 




16 


4 


16 


6 


128 


64 


96 


1 




17 


4 


16 


6 


32 


64 


96 


1 




18 


4 


16 


6 


64 


128 


96 


1 




19 


4 


16 


6 


64 


32 


96 


1 




20 


4 


16 


6 


64 


32 


192 


1 




21 


4 


16 


6 


64 


32 


48 


1 




22 


4 




6 


256 




256 


1 


axisymm. 


23 


4 




6 


512 




512 


1 


axisymm. 


24 


4 




6 


256 




256 


0.05 


axisymm. 


25 


4 




6 


512 




512 


0.1 


axisymm. 


26 


4 




8 


256 




256 


0.2 


axisymm. 


27 


4 




6 


512 




512 


0.5 


axisymm. 


28 


4 




6 


256 




256 


2 


axisymm. 


29 


4 




6 


256 




256 


8 


axisymm. 



6.3.1 Vertical box size 

One concern in formulating a numerical model such as this is 
whether an artificial aspect of the boundary conditions con- 
trols the outcome. In this respect, the azimuthal and radial 
boundary conditions arc somewhat less worrisome than the 
vertical boundary conditions; the former restrict the scale of 
structure that can develop to the size of the box, while the 
latter introduces a reflection of outgoing waves that would 
otherwise dissipate in the disc atmosphere (although there 
could be a sharp boundary between a hot disc atmosphere 
and the body of the disc that would also produce reflections, 
albeit from a free rather than fixed boundary). 

We have tested for variation in linear growth rate with 
the vertical size of box: compare, e.g., runs 2 and 3. Both 
runs give (llui5k,z/dt ^ 0.250. Thus for values of close 
to the fiducial value, QH, we observe no change in the nu- 
merically measured growth rate. For the non-linear devel- 
opment, we observe a trend in that the decay of epicyclic 
energy is slower in the larger boxes. Quantitatively, approx- 
imately 20% more of the initial epicyclic energy is left in run 
5 at t = 75n~^ than in run 4 at the same instant. The small 
size of this effect suggests that the boundary condition is 
not controlling the outcome. 



6.3.2 Type of vertical boundary condition 

To test the effect of varying the character of the vertical 
boundary conditions, we performed a series of axisymmetric 



runs with periodic vertical boundary conditions (runs 6 to 
11). These runs arc special in that the initial conditions have 

Vx = S{LzI2tt) sm{2'Kz/Lz) + 6 (65) 

to make them compatible with the boundary conditions. 
This velocity profile then requires larger than reflect- 
ing boundary condition runs so as to better approximate a 
linear shear near the mid-plane of the disc; this in turn im- 
plies lower numerical resolution Uz/Lz for a given grid size 

Hz. 

The growth rates measured in the periodic vertical runs 
are not distinguishably different from those measured in re- 
flecting boundary condition runs at similar resolution. The 
decay of epicyclic energy also follows a similar trajectory un- 
der the two different boundary conditions. This again sug- 
gests that the boundary conditions are not controlling the 
outcome. 



6.3.3 Azimuthal box size 

One might be concerned that the limited azimuthal size of 
the box is somehow limiting the linear development of the 
instability or the non-linear outcome. This concern turns out 
to be justified over a limited range in Ly for both the linear 
and non-linear development of the instability. 

Figure 4 shows the evolution of NAi.2,3 in the fidu- 
cial run 1. This demonstrates that m = 1 structure grows 
most quickly, and suggests that the low-m modes avail- 
able in a larger box might grow even more quickly. Fig- 
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Figure 5. Colour image of density on slices through the fiducial 
model. The top panel shows an a; — j; slice at 2 = 0, the middle 
panel shows an x — z slice at j; = 0, and the bottom panel shows a 
y — z slice at a; = 0. Black is highest density, while blue is lowest. 



ure 6 demonstrates that this is indeed the case. It shows 
the evolution of NAi/Ly in runs 12, 13, 14, and 15, which 
have Ly = 8, 16, 32, 64 respectively. The numerical resolu- 
tion Uy/Ly is constant between these runs. The early growth 
of NAi is strongest for Ly = 64. In fact, the growth rate 
of NAi for Ly = 64 is 0.3211, close to the 0.42$! measured 
for dlni?k,z/dt. This evidence suggests that there is a non- 
axisymmetric counterpart of the parametric instability for, 
and only for, sufficiently low m/Ly. 

In retrospect this result is easy to understand. Our lin- 
ear analysis could be trivially generalized to weakly non- 
axisymmetric disturbances using the shearing wave formal- 
ism of non-axisymmetric density wave theory (Goldreich & 
Lynden-Bell 1965). The radial wavenumber of the distur- 
bance would then evolve according to kx = qflky. The am- 
phtudes of these shearing waves will grow so long as the 
radial wavenumber does not change so much over a growth 
time that the resonance between the daughter modes and 
parent mode is detuned. 

The non-Unear development of the instability is not 
strongly altered by variation of Ly. There is a tendency for 
Ecpi to decay slightly more slowly in models with larger 
Ly. Presumably this is because large-scale, slowly dissipat- 




tn 

Figure 6. Evolution of NAi/L^ in runs with Ly = 8,16,32, 
and 64. The ordinate units are arbitrary. 

ing motions are available in the large Ly box, and these 
motions act as a reservoir for the epicyclic energy. 

6.3.4 Numerical resolution 

Runs 16 through 21, 1, 22, and 23 test the effect of numerical 
resolution. The measured linear growth rate is slightly de- 
pendent on resolution in the meridional plane. For nx,nz = 
32,48, 7num = dlnSk.^/dt = 0.37^; for n^^n^ = 64,96, 
7num = 0.460; for n^,n^ = 128,192, 7^^^ = 0.49!r2; for 
nx,nz — 512,512, 7num — 0.53f2. This should be compared 
to the linear theory rate 7 = 0.650. The non-linear decay 
rate of the warp also has a trend with resolution in that high 
resolution models decay more slowly. 

6.3.5 Initial warp amplitude 

At small S/Q., the three-mode coupling analysis of §5 pre- 
dicts that 

0.655. (66) 

To test this, we measured a numerical growth rate for a 
series of values of S. In each case, we used the run with 
the highest available numerical resolution in the meridional 
plane (we will be concerned only with the development of 
axisymmetric instabilities here). The results are listed in Ta- 
ble 2. The numerical growth rates typically lie ~ 15% below 
the analytic rates. We believe that this is primarily a resolu- 
tion effect, since the growth rate grows monotonically with 
increasing resolution. 

At large S/Q the three- mode coupling analysis is no 
longer applicable. The incompressible, unstratified model of 
§4 is not limited to small S/fl, and numerical integrations 
of those model equations show reduced growth for S/Q > 1. 
This is observed in the numerical models listed in Table 2. 
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Table 2. List of measured and predicted values for growth rate 



No. 


s/n 


O'num / ^ 


7th /f^ = 3V3S/S 


7num/7t 


24 


0.05 


0.020 


0.032 


0.62 


26 


0.1 


0.055 


0.065 


0.85 


26 


0.2 


0.11 


0.13 


0.85 


27 


0.5 


0.28 


0.32 


0.87 


23 


1.0 


0.49 


0.65 


0.76 


28 


2.0 


0.82 


1.3 


0.63 


29 


8.0 


1.7 


5.2 


0.33 



7 DISCUSSION 

We have shown that warps in unmagnetized, non-self- 
gravitating, inviscid Keplerian discs are subject to an in- 
stability that leads, in the non-linear regime, to damping 
of the warp. The damping is non-linear in that the growth 
rate of the instability is proportional to the amplitude of the 
warp. How might this process influence the development of 
warps in astrophysical discs? 

A full answer to this question awaits an understanding 
of other damping mechanisms that might act on a warp. 
If the disc is already turbulent, the parametric instability 
may not occur. In the Appendix, we estimate the influence 
of an isotropic effective turbulent viscosity on our analysis. 
We find that the instability still occurs if the amplitude of 
the in-plane oscillations is sufficiently large, possibly S/Q ^ 
30q, where a is the dimensionless viscosity parameter. In 
addition, for a given warping of the disc, is reduced in 
the presence of viscosity because the resonance that drives 
the in-plane oscillations is weakened. Therefore our analysis 
is most relevant to discs with a <^ 1. 

Assuming that the instability develops into the non- 
linear regime, we may hazard a guess at the consequences 
for the warp, using a simple model for the interaction of in- 
plane and out-of-plane motions in the warp. We then apply 
the result to the warped, masing disc in NGC 4258. 

Consider a linear, axisymmetric bending wave in the 
corotating frame. The radial and vertical displacements as- 
sociated with the wave, ~ 2:exp(ifer) and ~ exp(ifer), 
obey 



d^ 

dt2 



+ Q ^2 ^ ifcCs^r^zCr 



(67) 



(68) 



It is evident that the in-plane and out-of-plane motions 
behave as coupled harmonic oscillators. If k ~ and 
H/r <^ 1, as we shall assume, then the coupling is weak 
in the sense that kcs <C fi. The coupling would be negligible 
except that ~ il^ for a Keplerian disc, so the vertical and 
radial oscillations are resonant. 

We will use equations ( |67|6^ as the basis of a simple 
model for the warp. In addition to the terms from linear 
theory, we will include an out-of-plane driving term, meant 
to model the Pringle instability, and a non-linear damping 
term, meant to model the effects of the parametric instabil- 
ity. In terms of the normalized displacements = £,z/H 
and Xr = i^r{z = H)/H, the governing equations are 



-e^'^Xr^-lQ.^^ 



at' 



(70) 



Here e = kc^/Q., 7/2 is a dimensionless growth rate for the 
driving instability, and / is some function of Xr = dXr/&t 
that models the damping effect of the parametric instability 
on in-plane motions. 

We set n = Q, = 1 and look for steady-state oscilla- 
tory solutions of equations (^ To start, define "action- 
angle" variables by 



X, 

Xr 



a 2 cos 
dr cos 



flzCOS^z, Xz 
a,. COS 6r , Xr 

Consistency of these definitions requires 
6rar sin Or - 



-Uz sm f 
-ttr sin t 



~az sm Oz, 
~ar sin 6r. 



(71) 
(72) 

(73) 
(74) 



These relations and the differential equations for Xr and Xz 
yield 



702 sin 6z + ear cos Or sin 
/ sin Or + eaz sin Or cos Oz , 
1 -I- 7 sin Oz cc 



H cos t 

fflz 



1 -I- — cos Or H cos Or COS f 



(75) 
(76) 
(77) 

(78) 



When 7, e, and / are <^ 1, the rapidly-oscillating terms can 
be averaged out: 



A6' 



/Oz _ Or \ 

\ar az) 



cos AS, 



-ar sin A9, 



/(ar) + -ttz sin A6I, 



(79) 
(80) 
(81) 



where A6 = Or — Oz, f = (fsinOr), and {) denotes a phase 
average. Since / should be an odd function of Xr, we have 

if cos Or) = 0. 

A steady state requires either Uz = ar or cos AO = 0. In 
astrophysical applications the latter solutions are likely to be 
more relevant, because thin discs with observable warps have 
az ^ 1, while ar ^ 1 would require implausibly supersonic 
motions within the disc. 

The steady state solutions with cos A^ = have 



fjar) ^ 
ar 27 ' 

az = (e/7)ar. 

Stability of this solution requires 

ar < az 

and, for n = d In //d In a^, 
-hr <n<l. 



(82) 
(83) 

(84) 
(85) 



(69) 



To go further, we must guess at /(flr). At very small ar, 
we expect n > 2 because the growth rate of the paramet- 
ric instability is proportional to a^. At larger amplitudes / 
models the effect of stationary, fully developed turbulence. 
If turbulence initiated by the parametric instability acts as 
a viscosity (as might also be expected in an MHD turbulent 
disc) then n ~ 1. Notice that n ~ 1 is consistent with the 
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approximately exponential damping of the warp observed 
in the numerical experiments of §6. The decay time of the 
epicyclic energy shown in Figure 2 is ~ 30$!"^, which could 
be fitted by / ~ a^/GO: i.e., by a dimensionless viscosity 

Let us apply these results to the masing disc in NGC 
4258 (Miyoshi et al. 1995, also Gammie et al. 1999 and 
references therein). The presence of water vapour sug- 
gests Cs ~ Ikmsec"^. The rotational velocity is Vc ~ 
1100(r/0.13pc)~^''^kmsec^^ Thus e ~ H/r « 10"^ The 
orientation of the inner edge of the masing disc at ~ 0.13 pc 
differs from that at the outer edge at ~ 0.24 pc by ~ 0.25 
radians, so ~ 250. The corresponding velocity of in-plane 
motions at z — H would be ~ 250 km see" ^, were linear 
theory applicable. Finally, an estimate of the growth time of 
the Pringle instability gives r ~ 5 x 10® yr ( Maloney et al. 
199e|)7^r 7 = 2/(nT) ~ 1.2 x 10"* if Q. is taken at the outer 
edge of the disc. Direct measurements of maser velocities 
show that the departure from a Keplerian rotation curve is 
small, so we may take «; = fi. 

Of the model parameters e, flz, and 7, the last is most 
uncertain. Let us assume that the disc is described by our 
model equations and is in equilibrium, then solve for 7. 
Suppose / « asign{Xr)\Xr\"' / X^~^ , where Xq is a char- 
acteristic amplitude, and n is just slightly smaller than 1 
so the equilibrium is stable. Here a is the dimensionless 
strength of the damping. For n = 1 the damping rate is 
precisely what one would expect for a viscous disc with 
1/ — QfCg/n. Then / — — Qp„+isign(ar)|arl"/^o'~^, where 
Pn+i = (I sin6lr-r+^), and 



-l/n 1 + 1/n 

e Hz 



-1 + 1/n -vrl- 
^0 



-1/r. 



(86) 

This is smaller 



7 = (2ap„+i) 

In the limit ri ^ 1, 7 ^ c'^/a — 10~''/a. 
than the nominal value by a factor of 100a. 

We can obtain an upper limit on 7 if we suppose that 
the warp is limited by coupling to in-plane motions and 
that ar < 1. Then, to satisfy equation (^3|), we must have 
7 < e/oz « 4 X 10~®. This limit is due to the weakness 
of coupling between in-plane and out-of-plane motions. For 
larger growth rates the warp must be limited by some other 
process, or perhaps it is not limited and the disc is destroyed. 
The maximum dissipation rate to be expected of any local 
mechanism (one that relies upon the internal dynamics of 
the disc) is ~ c^n per unit mass. The local energy in verti- 
cal motions of the warp as measured in the corotating frame 
of the gas is E„ ~ per unit mass, and the rate at 

which the Pringle instability adds to this energy is 



It is therefore difficult to see how any local mechanism can 
limit the warp if the growth rate is as large as Maloney et 
al.'s estimate. 

Part of this discrepancy might be explained by uncer- 
tainties in the parameters of Maloney et al.'s model. It is 
also possible that non-linear effects reduce the growth rate 
of the Pringle instability. There is some evidence for this in 
that the disc in NGC 4258 is not as strongly twisted as one 
would expect from linear theory (see Herrnstein, 1997 for a 
warp model that fits the observations). Finally, non- linear 
effects may also enhance the coupling e between in-plane 
and out-of-plane oscillations. 
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APPENDIX A: EFFECT OF VISCOSITY ON 
THE INSTABILITY 

Here we make a rough estimate of the influence of pre- 
existing turbulence on the development of the parametric 
instability. In the absence of a proper understanding of the 
interaction of the in-plane oscillations with turbulence (see 
Torkelsson et al. 2000 for a direct numerical approach) we 
suppose that the process may be treated as an isotropic ef- 
fective viscosity of the form 



(Al) 



where a is the usual dimensionless viscosity parameter. 

The Floquet analysis of §4 is easily adapted to include 
viscosity. The solution is identical but multiplied by the 
damping factor 
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exp 



(A2) 



k, 30c 



(A5) 



When S/Q is small the time-dependence of k may be ne- 
glected and the growth rate of the instability is simply re- 
duced according to 

s 1-^ s — uk^ . (A3) 

For a tuned mode with k^/kz = ±\/3, the modified growth 
rate is 

s^^S-^{k.Hfa^. (A4) 

In the unstratified model k^ is unrestricted, but, in reality, 
the mode must fit inside the disc. An approximate criterion 
can be obtained by requiring that at least half a wavelength 
fit into a vertical distance 2H, i.e. k^H 7r/2. The instabil- 
ity then develops provided that 

S_ 

The efi'ect of a small viscosity on the mode-coupling 
analysis of §5 can also be estimated. Provided that bulk 
viscosity effects may be neglected, the eigenfunctions of wave 
modes in an isothermal disc (§3) can still be obtained in 
terms of Hermite polynomials when viscosity is included. 
The dispersion relation is much more complicated, but can 
be expanded for small a to obtain 

= J°\k) - iQM„(fc) -I- 0{a), (A6) 

where the first term is the frequency in the inviscid case 
(equation ^) and d„ is a dimensionless measure of the vis- 
cous damping rate of the mode. Generally speaking, d„ in- 
creases with increasing n and k. We therefore consider the 
first parametric resonance, which occurs between modes 1 
and 2 (i.e. tJi -I- aJ2 = f2) at kH ~ 1.88. The inviscid growth 
rate is then s ~ 0.1385. At this point di « d2 ~ 4.5, giving a 
net growth rate of 0.1385 — 4.50^. The instability proceeds 
only if S/^l k 33q, which agrees closely with equation (A5). 
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